module cmfdmg
    use mtblmg
    implicit none
    real(8), private, pointer, dimension(:,:,:)     :: phi
    real(8), private, pointer, dimension(:,:)       :: xfdmtot, xfdmrmv, xfdmnuf, xfdmchi
    real(8), private, pointer, dimension(:,:,:)     :: xfdmsct
!    real(8), private, pointer, dimension(:,:,:)     :: tot, rmv, nuf, dif
!    real(8), private, pointer, dimension(:,:,:,:)     :: sct
    real(8), private, pointer, dimension(:,:)       :: psi, psip
    type(mbmg), private :: mtcmfdmg
    real(8), private, pointer, dimension(:,:,:,:)   :: dtilde, dhat
    integer, private :: maxout, maxin 
    real(8), private :: ccout, ccin
    integer, private :: nx_, ny_, nxy_, nz_, ng_
    logical, private :: l1gcmfd, ifwldt
    real(8), private, pointer, dimension(:,:,:) :: phi_adj
    integer, private, pointer, dimension(:) :: iix, jjy ! indexing from center
    real(8), private, pointer, dimension(:,:,:)   :: q
    integer, private, pointer, dimension(:,:) :: map_asmcal
    
!    logical, private, pointer, dimension(:,:,:,:) :: 
contains

subroutine initcmfdmg
    use allocs
    use asmcal
    use bicgmg, only:mallocbicgmg
    include "pspec.FI"
    include "xsec.FI"
    include "hmgz.FI"
    integer :: xm, ym, i, j
    ifwldt = lwldt
	nx_=nfmx;ny_=nfmy;nxy_=nfmx*nfmy;nz_=1
    ng_=ng
    if (MC_ASMCAL) ng_=ngcd
    call dmalloc0(phi, 0,nxy_, 0,nz_, 1,ng_)
    call dmalloc(psi, nxy_, nz_)
    call dmalloc(psip, nxy_, nz_)
    call dmalloc(dtilde, 4, ng_, nx_, ny_)
    call dmalloc(dhat, 4, ng_, nx_, ny_)

    call dmalloc(phi_adj, 2, nxy_, nz_)
    call dmalloc(q, nxy_,nz_,ng_)

!    tot=>htot
!    rmv=>hrmv
!    nuf=>hnuf
!    sct=>hmsct
!    dif=>hdif
    call mallocbicgmg(nxy_, nx_, ny_, nz_, ng_)
    call newmbmg(mtcmfdmg, nxy_, nx_, ny_, nz_, ng_)
    maxout = max_out
    maxin  = max_in
    ccout = err_out
    ccin  = err_in
    
    call dmalloc(iix, nx_)
    call dmalloc(jjy, ny_)
    xm = int(nx_/2)
    ym = int(ny_/2)
    do i=1, xm;iix(i)=xm-i+1;enddo
    do i=xm+1, nx_; iix(i)=i;enddo
    do j=1, ym;jjy(j)=xm-j+1;enddo
    do j=ym+1, ny_; jjy(j)=j;enddo
    

end subroutine

subroutine initFDMMG
    use allocs
    include "pspec.FI"
    include "xsec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    integer :: i, j, g, l, tpin, ii, tc, gp, tpm
    real(8) :: rvp
    type (strcpin), pointer :: ppin
    
    
    rvp=1._8/vpin
    
    ! set uniform flux
    phi=1._8
    
    
    if (MC_ASMCAL) then
        call mapping_asmcal()
    else
        call dmalloc(xfdmtot, ng_,ntypepin)
        call dmalloc(xfdmrmv, ng_,ntypepin)
        call dmalloc(xfdmnuf, ng_,ntypepin)
        call dmalloc(xfdmsct, ng_,ng_,ntypepin)
        call dmalloc(xfdmchi, ng_,ntypepin)
        ! generate homogenized cross section
        do i=1, ntypepin
            ppin=>typepin(i)
            do g=1, ng_
                do l=0, ppin%nr
                    tc=ppin%rc(l)
                    xfdmtot(g,i)=xfdmtot(g,i)+ppin%rv(l)*xtot(tc,g)
                    xfdmnuf(g,i)=xfdmnuf(g,i)+ppin%rv(l)*xnuf(tc,g)
                    do gp=1, ng_
                        xfdmsct(gp,g,i)=xfdmsct(gp,g,i)+xsctm(gp,g,tc)*ppin%rv(l)
                    enddo 
                    if (l .eq. 0) then
                        xfdmchi(g,i)=xchi(tc,g)
                    endif
                enddo
                xfdmtot(g,i)=xfdmtot(g,i)*rvp
                xfdmnuf(g,i)=xfdmnuf(g,i)*rvp
                do gp=1, ng_
                    xfdmsct(gp,g,i)=xfdmsct(gp,g,i)*rvp
                enddo 
                xfdmrmv(g,i)=xfdmtot(g,i)-xfdmsct(g,g,i)
            enddo
        enddo
    endif
!    do g=1, ng_
!        chi(g)=xchi(1,g)
!    enddo
    

!    do j=1,ny_;do i=1,nx_;
!        tpm=cplst(i,j)
!        if (tpm.eq.0) cycle ! for vacuum assembly
!        htot(:,i,j)=xfdmtot(:,tpm)
!        hrmv(:,i,j)=xfdmrmv(:,tpm)
!        hnuf(:,i,j)=xfdmnuf(:,tpm)
!        hmsct(:,:,i,j)=xfdmsct(:,:,tpm)
!        hdif(:,i,j)=1._8/3/xfdmtot(:,tpm)
!    enddo;enddo
contains
subroutine mapping_asmcal
    use asmcal 
    include "pspec.FI"
    integer:: ja, jp, ia, ip, tc
    integer:: i, j
    allocate(map_asmcal(nfmx, nfmy))
    map_asmcal=0
    
    j=0
    do ja=1, nay
    do jp=1, npy
        j=j+1
        i=0
        do ia=1, nax
        do ip=1, npx
            i=i+1
            tc=xsec_map(ip,jp, core(ia,ja))
            map_asmcal(i,j)=tc
        enddo
        enddo
    enddo
    enddo

endsubroutine

endsubroutine

subroutine makeMT(ldhat)
    use asmcal
    use define
    include "pspec.FI"
    include "xsec.FI"
    integer :: i, j, k, l, g, ia, ja, dr, lo, oi, oj, tc
    integer :: nxs, nxe
    real(8) :: hfm(4), dt, dh
    logical :: ldhat
    
    hfm(1:2) = hfmx
    hfm(3:4) = hfmy
    
    do k=1, 1
    do g=1, ng_; l=0
    do j=1, ny_;ja=f2cy(j);
        nxs=0;nxe=1
        do i=1, nx_;ia=f2cx(i);
            l=l+1
            if (MC_ASMCAL) then
                tc=map_asmcal(i,j)
            else
                tc=cplst(i,j)
            endif
            if (tc .eq. 0) cycle
            if (0 == nxs) nxs=i
            nxe=i
            
            call setnodel(mtcmfdmg, i,j,l)
            if (MC_ASMCAL) then
                call setmbmg(mtcmfdmg,  l,k,0,g, xsec_set.rmv(g,tc)*vfm)
            else
                call setmbmg(mtcmfdmg,  l,k,0,g, xfdmrmv(g,tc)*vfm)
            endif
            if (ldhat) then ! with dhat
                do dr=1, 4
                    lo=neighf(i,j,oi,oj,dr)
                    if (lo .eq. 0) then
!                        dt=dtilde(dr,g,i,j)
                        dh=dhat(dr,g,i,j)
                        call accmbmg(mtcmfdmg,l,k,0, g,    dh*hfm(dr)   )
                    else
                        dt=dtilde(dr,g,i,j)
                        dh=dhat(dr,g,i,j)
                        call accmbmg(mtcmfdmg,l,k,0, g,    (dt-dh)*hfm(dr)   )
                        call setmbmg(mtcmfdmg,l,k,dr,g,   -(dt+dh)*hfm(dr),lo)
                    endif
                enddo  
            else
                do dr=1, 4
                    lo=neighf(i,j,oi,oj,dr)
                    if (lo .eq. 0) then
                        dt=dtilde(dr,g,i,j)
                        call accmbmg(mtcmfdmg,l,k,0, g,    dt*hfm(dr)   )
                    else
                        dt=dtilde(dr,g,i,j)
                        call accmbmg(mtcmfdmg,l,k,0, g,    dt*hfm(dr)   )
                        call setmbmg(mtcmfdmg,l,k,dr,g,   -dt*hfm(dr),lo)
                    endif
                enddo  
            endif
        enddo
        call setnxse(mtcmfdmg,j,nxs,nxe) 
    enddo
    enddo 
    enddo
    
    call facilumg(mtcmfdmg)

end subroutine




subroutine calDtilde ! ==========================================================================
    use define
    use asmcal  
    include "pspec.FI"
    include "xsec.FI"
    include "hmgz.FI"
    integer :: i, j, k, l, g, ia, ja, dr, lo, oi, oj, tc
    integer :: nxs, nxe
    real(8) :: hfm(4), dt
    
    hfm(1:2) = hfmx
    hfm(3:4) = hfmy
    
    do k=1, 1
    do g=1, ng_; l=0
    do j=1, ny_
        do i=1, nx_
            l=l+1
            if (MC_ASMCAL) then
                tc=map_asmcal(i,j)
            else
                tc=cplst(i,j)
            endif
            if (tc .eq. 0) cycle
            do dr=1, 4
                lo=neighf(i,j,oi,oj,dr)
                if (lo .eq. 0) then
                    dtilde(dr,g,i,j)=dtildebndry(i,j,g,dr)
                else
                    dtilde(dr,g,i,j)=dtildeinterface(i,j,oi,oj,g)
                endif
            enddo  
        enddo
    enddo
    enddo 
    enddo
contains
function dtildebndry(mi,mj,g, dr)
    integer, intent(in) :: mi,mj,g, dr
    real(8) :: dtildebndry, totm, beta
    if (MC_ASMCAL) then
        totm=xsec_set.tot(g,map_asmcal(mi,mj))
    else
        totm=xfdmtot(g,cplst(mi,mj))
    endif
    beta=1._8/3/totm/hfmx
    dtildebndry= 2*albedo(dr)*beta/(albedo(dr)+2*beta)    
endfunction

function dtildeinterface(mi,mj,oi,oj,g)
    integer, intent(in) :: mi,mj,oi,oj,g
    real(8) :: dtildeinterface, totm, toto
    ! assuming that hfmx=hfmy
    if (MC_ASMCAL) then
        totm=xsec_set.tot(g,map_asmcal(mi,mj))
        toto=xsec_set.tot(g,map_asmcal(oi,oj))
    else
        totm=xfdmtot(g,cplst(mi,mj))
        toto=xfdmtot(g,cplst(oi,oj))
    endif
    dtildeinterface =2._8/(totm+toto)/3/hfmx
end function 
    
endsubroutine ! ==========================================================================

!subroutine calDhat
!    use define
!    include "pspec.FI"
!    include "xsec.FI"
!    include "hmgz.FI"
!    integer :: i, j, k, l, g, ia, ja, dr, lo, oi, oj, tc
!    integer :: nxs, nxe
!    real(8) :: hfm(4), dt
!    
!    hfm(1:2) = hfmx
!    hfm(3:4) = hfmy
!    
!    do k=1, 1
!    do g=1, ng; l=0
!    do j=1, ny_
!        do i=1, nx_
!            tc=cplst(i,j)
!            if (tc .eq. 0) cycle
!            do dr=1, 4
!!                if (hj(dr,g,i,j) .eq. 0._8) cycle
!                lo=neighf(i,j,oi,oj,dr)
!                if (lo .eq. 0) then
!                    dhat(dr,g,i,j)=dhatbndry(i,j,g,dr)
!                else
!                    dhat(dr,g,i,j)=dhatinterface(i,j,oi,oj,g,dr)
!                endif
!                
!            enddo  
!        enddo
!    enddo
!    enddo 
!    enddo
!contains
!function dhatbndry(mi,mj,g, dr)
!    include "hmgz.FI"
!    integer, intent(in) :: mi,mj,g, dr
!    real(8) :: dhatbndry
!    dhatbndry =0._8
!    if(hphi(g,mi,mj) .ne. 0._8) then
!        dhatbndry = hj(dr,g,mi,mj)/hphi(g,mi,mj)
!    else
!        Dhatbndry = Dtilde(dr,g,mi,mj)
!    endif
!endfunction
!
!function dhatinterface(mi,mj,oi,oj,g, dr)
!    include "hmgz.FI"
!    integer, intent(in) :: mi,mj,oi,oj,g, dr
!    real(8) :: dhatinterface
!    real(8) :: netj
!    netj = hj(dr,g,mi,mj)
!    if (netj .ne. 0._8) then
!        dhatinterface = -( netj + Dtilde(dr,g,mi,mj)*(hphi(g,oi,oj)-hphi(g,mi,mj)) )/(hphi(g,mi,mj)+hphi(g,oi,oj))
!    else
!        dhatinterface = 0._8
!    endif
!end function 
!    
!endsubroutine ! ==========================================================================

subroutine readfdm()
	open (7001, file='fdm.bin', form='binary', status='old')
	read (7001) psi
	close(7001)
	call genFSPDF
endsubroutine

subroutine genFSPDF
    include "arrays.FI"
    include "pspec.FI"
    integer :: l, i, j
    real(8) :: sumadj, rmsfdm, relerr(nx_)
    l=0
	sumadj=0
    do j=1,ny_;do i=1, nx_
        l = l+1
        if (cplst(i,j) .eq. 0) then ! for vacuum assembly
            initFSD(i,j) = 0._8
        else
            initFSD(i,j)=psi(l,1)
            sumadj = sumadj + initFSD(i,j)
        endif
    enddo;enddo
    do j=1,ny_;do i=1, nx_
        initFSD(i,j) = initFSD(i,j)/sumadj
    enddo;enddo

    open (7001, file='fdmerr.txt', status='unknown')

    ! evaluate the error of given distribution
    if (ireft .eq. 1) then 
        rmsfdm=0
        do j=1, nfmy
            relerr=0
            do i=1, nfmx
                if (refpsi(i,j)>0) relerr(i)=(initFSD(i,j)-refpsi(i,j))/refpsi(i,j)
                rmsfdm=rmsfdm+(initFSD(i,j)-refpsi(i,j))**2
            enddo
            write(7001,'(2000(es12.5,1x))') relerr
        enddo
    endif
	close(7001)
    rmsfdm = sqrt(rmsfdm*npowpinref)
    write (*, '("RMS Error of FDM FSD :", es12.5)'), rmsfdm
    
endsubroutine

!subroutine cmfdmg_mc(k, err, iter, err_ref)
!    use define
!    include "pspec.FI"
!    include "hmgz.FI"
!    integer, intent(out)   :: iter
!    real(8), intent(inout) :: k
!    real(8), intent(out)   :: err, err_ref
!    real(8) :: dr, nrm, sumadj
!    logical :: ldhat
!    integer :: i, j, l, g, nzero, ii, jj
!    integer :: nnz, lo, oi, oj, dir, gp
!    ldhat = .true.
!    
!    ! copy homogenized flux 
!    
!    nzero=0
!    do g=1, ng_;l=0
!    do jj=1, ny_;j=jjy(jj)
!    do ii=1, nx_;i=iix(ii)
!        l=l+1
!        if (cplst(i,j) .eq. 0) cycle
!        if (hphi(g,i,j) > 0._8) then
!        else
!            ! FIXING HOMOGENIZED CROSS SECTIONS 
!            nnz = 0;
!            do dir=1, 4
!                lo=neighf(i,j,oi,oj,dir)
!                if (lo .eq. 0) cycle
!                if (hphi(g,oi,oj)>0) then
!                    nnz=nnz+1
!                    hphi(g,i,j)=hphi(g,i,j)+hphi(g,oi,oj)
!                    htot(g,i,j)=htot(g,i,j)+htot(g,oi,oj)
!                    hnuf(g,i,j)=hnuf(g,i,j)+hnuf(g,oi,oj)
!                    hphi(g,i,j)=hphi(g,i,j)+hphi(g,oi,oj)
!                    hRmv(g,i,j)=hRmv(g,i,j)+hRmv(g,oi,oj)
!                    do gp=1, ng_
!                        hmsct(gp,g,i,j)=hmsct(gp,g,i,j)+hmsct(gp,g,oi,oj)
!                    enddo
!                endif
!            enddo
!            if (nnz .eq. 0) then ! severe error
!!                print*, "no infromation from other nodes"
!!                cycle
!            nzero = nzero+1
!            else
!                hphi(g,i,j)=hphi(g,i,j)/nnz 
!                htot(g,i,j)=htot(g,i,j)/nnz 
!                hnuf(g,i,j)=hnuf(g,i,j)/nnz 
!                hphi(g,i,j)=hphi(g,i,j)/nnz 
!                hRmv(g,i,j)=hRmv(g,i,j)/nnz 
!                hmsct(:,g,i,j)=hmsct(:,g,i,j)/nnz
!            endif
!            
!!            do j=1, ncmy;do i=1, ncmx;do g=1, ngcmfd;do k=1, 4
!!                hj(k,g,i,j)=hjac(k,g,i,j)-hjac(LTOPST(k),g,i+LTNEIGH(1,k),j+LTNEIGH(2,k))
!!            enddo;enddo;enddo;enddo
!            
!        endif
!        
!        phi(l,1,g)=hphi(g,i,j)
!    enddo
!    enddo
!    enddo
!    
!    if(nzero>0) then
!        print *, "Unsafe to perform MG-CMFD. There are (", nzero, ") untracked meshes."
!        return
!    endif
!
!    call calDtilde
!    call calDhat
!    call makeMT(ldhat)
!    call poweriter(ldhat, k, iter, err, dr, .true., .false.)
!    
!    call genFSPDF
!    
!    if (ireft .eq. 1) then
!        nrm = 1._8/sum(psi)
!        err_ref=0._8
!        do j=1, ny_;do i=1, nx_
!            err_ref=err_ref+(refpsi(i,j)-psi(i,j)*nrm)**2
!        enddo;enddo
!        err_ref=sqrt(err_ref*npowpinref);
!    endif
!    
!        ! check the solution
!    do g=1, ng_
!    do l=1, nxy_
!        if ( psi(l,1)<0) then
!            print *, "negative flux"
!            iter = 0
!            err  = 1.0_8
!            do i=1, nxy_
!                phi(:,1,g) = 1._8
!            enddo
!            fdbsafe = .false.
!            return
!        endif
!    enddo
!    enddo
!    fdbsafe = .true.
!!     hphi = 0
!
!endsubroutine
!

subroutine FDM
    use cmfd2g, only :cmfd2g4mg
    include "pspec.FI"
    integer :: iter, i, j, g
    real(8) :: k, err, dr, sumadj
    logical :: lflog
    integer :: iter_cmfd
    real(8) :: err_cmfd
    real(8) :: sumpsi, rnrm
    integer :: npsi, l, st, ed
    lflog=.true.
    k=1._8
    call calDtilde
    do g=1,ng_;do i=1, nxy_
        phi(i,1,g)=1._8
    enddo;enddo
!    err = err_out
!	iter = max_out
    err = 1e-3_8
	iter = 0
    call makeMT(.false.)
    call poweriter(.false., k, iter, err, dr, .false., lflog)
!    call poweriter(.false., k, iter, err, dr, ifwldt, lflog)
!    call poweriter(k, iter, err, dr, .false., lflog)
    
    if (dopprec_fdm) call pinpower_reconstruct
    
    call genFSPDF
    
    sumpsi =0._8
    l=0;npsi=0
    do j=1, ny_;do i=1, nx_
        l=l+1
        sumpsi =sumpsi+psi(l,1)
        if (psi(l,1)>0) npsi=npsi+1
    enddo;enddo
    rnrm = 1._8/sumpsi*npsi
    
    l=0
    do j=1, ny_;do i=1, nx_
        l=l+1
        psi(l,1)=psi(l,1)*rnrm
    enddo;enddo
    
    open (7001, file="FDM.txt", status='unknown')
    do j=ny_, 1
        st = (ny_-1)*nx_+1
        ed = ny_*nx_
        write(7001, '(2000(es12.5,1x))'), psi(st:ed,1)
    enddo
    close(7001)
    open (7001, file='fdm.bin', form='binary', status='unknown')
    write (7001) psi
	close(7001)

!    open (7001, file='fdm.txt', status='unknown')
!    do j=1,ny_
!        write (7001,'(2000(es12.5,1x))') psi((j-1)*nx_+1:j*nx_,1)
!    enddo
!	close(7001)
contains

subroutine pinpower_reconstruct
    use asmcal, only : pprec
    integer :: ia, ja, jp
    integer :: xbase, ybase, base1d
    real(8), dimension(npx, npy) :: psipart
    do ja=1, nay
        do ia=1, nax
            psipart=0
            do jp=1, npy
                ybase=(ja-1)*npy+jp-1
                xbase=(ia-1)*npx
                base1d=ybase*nfmx+xbase
                psipart(:,jp)=psi(base1d+1:base1d+npx,1)
            enddo
            call pprec(psipart,ia,ja)
            do jp=1, npy
                ybase=(ja-1)*npy+jp-1
                xbase=(ia-1)*npx
                base1d=ybase*nfmx+xbase
                psi(base1d+1:base1d+npx,1)=psipart(:,jp)
            enddo
        enddo
    enddo
endsubroutine	
    
endsubroutine

subroutine CMFDMGFS 
    use asmcal
    include "pspec.FI"
    integer :: i, j, l, g, tpm, ia, ja, tc
    psi=0
    do g=1, ng_;l=0
    do j=1, ny_;ja=f2cy(j);
    do i=1, nx_;ia=f2cx(i);
        l=l+1
        if (MC_ASMCAL) then
            tc=map_asmcal(i,j)
            if (tc .eq. 0) cycle ! for vacuum assembly
            psi(l,1)=psi(l,1)+xsec_set.nuf(g,tc)*vfm*phi(l,1,g)
        else
            tc=cplst(i,j)
            if (tc .eq. 0) cycle ! for vacuum assembly
            psi(l,1)=psi(l,1)+xfdmrmv(g,tc)*vfm*phi(l,1,g)
        endif
        
            
!        tc=cplst(i,j)
!        if (tc .eq. 0) cycle 
!        psi(l,1)=psi(l,1)+xfdmnuf(g,tc)*phi(l,1,g)*vfm
    enddo 
    enddo
    enddo
end subroutine

subroutine CMFDMGQ(k,g,q)
    use asmcal
    include "pspec.FI"
    real(8), intent(in) :: k
    integer, intent(in) :: g
    real(8), dimension(:,:,:), intent(out) :: q
    integer :: i, j, l, gp, tpm, ia, ja, tc
    l=0
    do j=1, ny_;
    do i=1, nx_;
        l=l+1
        if (MC_ASMCAL) then
            tc=map_asmcal(i,j)
            if (tc .eq. 0) cycle ! for vacuum assembly
            q(l,1,g)=xsec_set.chi(g,tc)*psi(l,1)/k
            do gp=1, ng_
                if (gp == g) cycle
                q(l,1,g)=q(l,1,g)+xsec_set.sct(g,gp,tc)*phi(l,1,gp)*vfm
            enddo
        else
            tc=cplst(i,j)
            if (tc .eq. 0) cycle ! for vacuum assembly
            q(l,1,g)=xfdmchi(g,tc)*psi(l,1)/k
            do gp=1, ng_
                if (gp == g) cycle
                q(l,1,g)=q(l,1,g)+xfdmsct(g,gp,tc)*phi(l,1,gp)*vfm
            enddo
        endif
        

!        tc=cplst(i,j)
!        if (tc .eq. 0) cycle ! for vacuum assembly
!        q(l,1,g)=xfdmchi(g,tc)*psi(l,1)/k
!        do gp=1, ng_
!            if (gp == g) cycle
!            q(l,1,g)=q(l,1,g)+xfdmsct(g,gp,tc)*phi(l,1,gp)*vfm
!        enddo
    enddo 
    enddo
end subroutine

subroutine makecurrent(ldhat) !===========================================================
    use define
    include "pspec.FI"
    include "xsec.FI"
    include "hmgz.FI"
    integer :: i, j, k, g, g2, dr, lo, oi, oj, l
    real(8) :: sj, dt, dh, po, pm
    logical :: ldhat
    
    hj2=0
    do g=1, ng_;g2=gc2(g);l=0
    do k=1, nz_
    do j=1, ny_
        do i=1, nx_
            l=l+1
            if (cplst(i,j) .eq. 0) cycle
            pm=phi(l,k,g)
!            if (ldhat) then
!                do dr=1, 4
!                    dh=dhat(dr,g,i,j)
!                    lo=neighf(i,j,oi,oj,dr)
!                    if (lo .eq. 0) then
!                        sj=dh*pm
!                    else
!                        dt=dtilde(dr,g,i,j)
!                        po=phi(lo,k,g)
!                        sj=-dt*(po-pm)-dh*(po+pm)
!                    endif
!                    hj2(dr,g2,i,j)=hj2(dr,g2,i,j)+sj
!                enddo
!             else
                do dr=1, 4
                    dt=dtilde(dr,g,i,j)
                    lo=neighf(i,j,oi,oj,dr)
                    if (lo .eq. 0) then
                        sj=dt*pm
                    else
                        po=phi(lo,k,g)
                        sj=-dt*(po-pm)
                    endif
                    hj2(dr,g2,i,j)=hj2(dr,g2,i,j)+sj
                enddo  
!            endif
        enddo
    enddo
    enddo
    enddo 
contains 
function current_bndry(mi,mj,g, dr)
    integer, intent(in) :: mi,mj,g, dr
    real(8) :: current_bndry, totm, beta
    totm=xfdmtot(g,cplst(mi,mj))
    beta=1._8/3/totm/hfmx
    current_bndry= 2*albedo(dr)*beta/(albedo(dr)+2*beta) *phi(l,1,g)   
endfunction
function current_interface(mi,mj,oi,oj,g)
    integer, intent(in) :: mi,mj,oi,oj,g
    real(8) :: dtilde, totm, toto, current_interface
    integer :: tpm, tpo
    ! assuming that hfmx=hfmy
    totm=xfdmtot(g,cplst(mi,mj))
    toto=xfdmtot(g,cplst(oi,oj))
    dtilde =2._8/(totm+toto)/3/hfmx
    current_interface = -dtilde*(phi(lo,1,g)-phi(l,1,g))
endfunction
endsubroutine makecurrent !===========================================================

subroutine Gcondense
    include "pspec.FI"
    include "xsec.FI"
    include "hmgz.FI"
    integer :: i, j, g, g2, gp, gp2, dr, k, l, tc
    real(8) :: phiseg, phisegr
    hphi2=0;
    htot2=0;
    hrmv2=0;
    hdif2=0;
    hnuf2=0;
    hmsct2=0;
    
    
    do g=1, ng_
    do k=1, 1
        g2=gc2(g)
        l=0
        do j=1, ny_
            do i=1, nx_
                l=l+1
                tc=cplst(i,j)
                if (cplst(i,j) .eq. 0) cycle
                ! group constants
                phiseg=phi(l,k,g)
                hphi2(g2,i,j)=hphi2(g2,i,j)+phiseg
                htot2(g2,i,j)=htot2(g2,i,j)+xfdmtot(g,tc)*phiseg
                hnuf2(g2,i,j)=hnuf2(g2,i,j)+xfdmnuf(g,tc)*phiseg
                do gp=1, ng_
                    gp2=gc2(gp)
                    hmsct2(gp2,g2,i,j)=hmsct2(gp2,g2,i,j)+xfdmsct(gp,g,tc)*phiseg
                enddo
            enddo
        enddo
    enddo
    enddo
    
    do g2=1, ng2
        do j=1, ny_
            do i=1, nx_
                if (cplst(i,j)  .eq. 0) cycle
                phisegr=1._8/hphi2(g2,i,j)
                htot2(g2,i,j)=htot2(g2,i,j)*phisegr
                hnuf2(g2,i,j)=hnuf2(g2,i,j)*phisegr
                do gp2=1, ng2
                    hmsct2(gp2,g2,i,j)=hmsct2(gp2,g2,i,j)*phisegr
                enddo
                hdif2(g2,i,j)=1._8/3/htot2(g2,i,j)
                hrmv2(g2,i,j)=htot2(g2,i,j)-hmsct2(g2,g2,i,j)
            enddo
        enddo
    enddo
endsubroutine

subroutine poweriter(ldhat, k, iter, err, dr, ifwldt, iflog)
    use cmfd2g, only : cmfd2g4mg, copypsi, cmfd2gacc
    use bicgmg, only : initbicg1g, solbicg1g
    use ifport
!$  use omp_lib
    implicit none
    include "xsec.FI"

    real(8), intent(inout)          :: k
    real(8), intent(inout)            :: dr, err
    integer, intent(inout)            :: iter
    integer                         :: i, j, g, gp, maxiterin, max_iter, iti, tot_in, l, gc
    real(8)                         :: errconvin, torErr, upper, upper2, lower, err_prv, ke, kd, kr
    real(8)                         :: r2, r20
    logical                         :: ldhat, ifwldt, iflog
    real(8)                         :: k_cmfd, err_cmfd, aa
	real(8)							:: errconv
	integer							:: iterconv
    integer                         :: iter_cmfd
                       
!$  real(8) :: ts(3)
!$  ts(3) = 0
!$      ts(1) = omp_get_wtime()
	print *, err, iter
	if (err .ne. 0._8) then
		errconv = err
	else
		errconv = ccout
	endif
	if (iter .ne. 0) then
		iterconv = iter
	else
		iterconv = maxout
	endif
    err    = 1.
	iter   =0
    dr = 1.
    if (ifwldt) kr=k
    tot_in=0                                

!    call makecurrent
!    call Gcondense
!    k_cmfd=k
!    call cmfd2gacc(k_cmfd, err_cmfd, iter_cmfd, psidiff)

    ! make fission source
    call cmfdmgFS

    do while (iter .lt. iterconv .and. err.gt.errconv)
        err_prv = err
        iter = iter + 1
        do g=1, ng_
            ! make total source
            call cmfdmgQ(k,g,q)
            call initbicg1g(mtcmfdmg,g,phi,q,r20)
            if (r20<1e-10) exit
            r2=r20;iti=0;
            do while (iti .lt. maxin .and. r2/r20.gt.ccin)
            !do while (iti .lt. 20)
!                print *, iti
                call solbicg1g(mtcmfdmg,g,r20,phi,r2)
                iti=iti+1
            enddo
!			print *, g, iti, r2, r20
20  format("    Inner Iteration | Group (", i3, ") completed.")
            if (g-g/10*10 == 0) print 20, g
        enddo
		if (ng_-ng_/10*10 .ne. 0) print 20, ng_
        tot_in = tot_in+iti
            
        ! make fission source from updated flux
        do i=1, nxy_
            psip(i,1) = psi(i,1)
        enddo
        call cmfdmgFS
                
        ! evaluate k                                            
        upper = 0.
        lower = 0.
        upper2= 0.
!!$omp parallel do private(i), reduction(+:upper,lower,upper2)       
        do i=1, nxy_
            upper = upper + psi(i,1)*psi(i,1)
            lower = lower + psi(i,1)*psip(i,1)
            upper2 = upper2 + (psi(i,1)- psip(i,1))*(psi(i,1)- psip(i,1))
        end do
!!$omp end parallel do

        k=k*upper/lower
                
        ! evaluate err 
        lower = sqrt(upper)
        upper = sqrt(upper2)
        err = upper/lower
		dr = err/err_prv  

10  format ("Outer Iteration Step(" i4, ") | k = ", f7.5, ",  err = ", es11.5, ",  dr = ", f7.5)
        print 10, iter, k, err, dr    
        

        ! two group CMFD with fine mesh
!        if(.false.) then
        if(ifwldt) then
            call makecurrent(ldhat)
            call gcondense
            k_cmfd=k
            iter_cmfd = 0
            err_cmfd  = 0
            call cmfd2g4mg(k_cmfd, iter_cmfd, err_cmfd, phi_adj)
            ! prolongation
            do g=1, ng_; gc=gc2(g)
                do l=1, nxy_
                    phi(l,1,g)=phi(l,1,g)*phi_adj(gc,l,1)
                enddo
            enddo
            !call cmfd2gacc(k_cmfd, err_cmfd, iter_cmfd, aa)
            k=k_cmfd 
!            call copypsi(psi)
        endif        
        
    end do
!$      ts(2) = omp_get_wtime()
    dr = err/err_prv  
    if(iflog) then
        write(*,*) iter, err, tot_in
        write(*, '(" Solution from MG CMFD : k=",f7.5,1x,"dr=",f7.5)'), k, dr
        write(*, '(" Computing time    : ", f12.3,1x,"(sec)")') ts(2)-ts(1)
    endif

!            call makecurrent(ldhat)
!            call gcondense
    
end subroutine


endmodule
